---
title: "conjoint analysis Yemen"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc)
require(stringi)
require(ggpubr)

dat <- read_csv("data/cj_data_Yemen.csv")
dat <- na.exclude(dat)

str(dat)
dat[,c(43:51)] <- data.frame(lapply(dat[,c(43:51)],as.factor))
dat[,c(53:55)] <- data.frame(lapply(dat[,c(53:55)],as.factor))

# filer finished sample
#dat <- dat |> 
#  filter(Finished == "1")

# Ceasefire: unexpected event during survey
dat$RecordedDate <- as.Date(stri_match_first_regex(dat$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat$ceasefire <- NA
dat$ceasefire[dat$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat$ceasefire[as.Date("2025-01-17") <= dat$RecordedDate] <- "Post"
dat$ceasefire <- factor(dat$ceasefire)

# task
dat$task <- as.factor(dat$task)

# Gender
dat <- dat |> 
  mutate(Gender = case_when(Gender == 1 ~ "Male",
                            Gender == 2 ~ "Female"))
dat$Gender <- as.factor(dat$Gender)

# Age group
dat <- dat |> 
  mutate(Age_group = case_when(D2 == 1 ~ "18-20",
                               D2 == 2 ~ "21-30",
                               D2 == 3 ~ "31-40",
                               D2 == 4 ~ "41-50",
                               D2 == 5 ~ "51-60"))
dat$Age_group <- as.factor(dat$Age_group)

# Education
dat <- dat |>
  mutate(Education = D3)
dat$Education <- as.factor(dat$Education)

# Income
dat <- dat |>
  mutate(Income = D5)
dat$Income <- as.factor(dat$Income)

# Sect
dat <- dat |> 
  mutate(Sect = case_when(D6 == 1 ~ "Sunni",
                          D6 == 2 ~ "Zaydi",
                          D6 == 3 ~ "Others"))
dat$Sect <- as.factor(dat$Sect)

# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q1_1 >= 67 ~ "High",
                              Q1_1 <= 66 & Q1_1 >= 34 ~ "Middle",
                              Q1_1 <= 33 ~ "Low"))
dat$Islamism <- as.factor(dat$Islamism)

# Secularism
dat <- dat |> 
  mutate(Secularism = case_when(Q1_2 >= 67 ~ "High",
                                Q1_2 <= 66 & Q1_2 >= 34 ~ "Middle",
                                Q1_2 <= 33 ~ "Low"))
dat$Secularism <- as.factor(dat$Secularism)

# Tribalism
dat <- dat |> 
  mutate(Tribalism = case_when(Q1_3 >= 67 ~ "High",
                               Q1_3 <= 66 & Q1_3 >= 34 ~ "Middle",
                               Q1_3 <= 33 ~ "Low"))
dat$Tribalism <- as.factor(dat$Tribalism)

# Arabism
dat <- dat |> 
  mutate(Arabism = case_when(Q1_4 >= 67 ~ "High",
                             Q1_4 <= 66 & Q1_4 >= 34 ~ "Middle",
                             Q1_4 <= 33 ~ "Low"))
dat$Arabism <- as.factor(dat$Arabism)

# US
dat <- dat |> 
  mutate(US = case_when(Q2_1 >= 4 ~ "Like",
                        Q2_1 == 3 ~ "Neither",
                        Q2_1 <= 2 ~ "Dislike"))
dat$US <- as.factor(dat$US)

# Iran
dat <- dat |> 
  mutate(Iran = case_when(Q2_2 >= 4 ~ "Like",
                          Q2_2 == 3 ~ "Neither",
                          Q2_2 <= 2 ~ "Dislike"))
dat$Iran <- as.factor(dat$Iran)

# Saudi
dat <- dat |> 
  mutate(Saudi = case_when(Q2_3 >= 4 ~ "Like",
                           Q2_3 == 3 ~ "Neither",
                           Q2_3 <= 2 ~ "Dislike"))
dat$Saudi <- as.factor(dat$Saudi)

# Turkey
dat <- dat |> 
  mutate(Turkey = case_when(Q2_4 >= 4 ~ "Like",
                            Q2_4 == 3 ~ "Neither",
                            Q2_4 <= 2 ~ "Dislike"))
dat$Turkey <- as.factor(dat$Turkey)

# UAE
dat <- dat |> 
  mutate(UAE = case_when(Q2_5 >= 4 ~ "Like",
                         Q2_5 == 3 ~ "Neither",
                         Q2_5 <= 2 ~ "Dislike"))
dat$UAE <- as.factor(dat$UAE)

# Israel
dat <- dat |> 
  mutate(Israel = case_when(Q2_6 >= 4 ~ "Like",
                            Q2_6 == 3 ~ "Neither",
                            Q2_6 <= 2 ~ "Dislike"))
dat$Israel <- as.factor(dat$Israel)

# Axis
dat <- dat |> 
  mutate(Axis = case_when(Q2_7 >= 4 ~ "Like",
                          Q2_7 == 3 ~ "Neither",
                          Q2_7 <= 2 ~ "Dislike"))
dat$Axis <- as.factor(dat$Axis)

# UN
dat <- dat |> 
  mutate(UN = case_when(Q2_8 >= 4 ~ "Like",
                        Q2_8 == 3 ~ "Neither",
                        Q2_8 <= 2 ~ "Dislike"))
dat$UN <- as.factor(dat$UN)


# Houthi Allied parties supporter
dat <- dat |> 
  mutate(Houthiallied_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Houthiallied_party_supporter = ifelse(Houthiallied_party_supporter >= 7, "Houthis Allied Party Supporter", "Others"))
dat$Houthiallied_party_supporter <- as.factor(dat$Houthiallied_party_supporter)

# Patronage
dat <- dat |> 
  mutate(Patronage = case_when(Q4 >= 3 ~ "High",
                               Q4 <= 2 ~ "Low"))
dat$Patronage <- as.factor(dat$Patronage)

# High patronage and Houthi allied supporter
dat <- dat |>
  mutate(Patron_Houthi = ifelse(Houthiallied_party_supporter == "Houthis Allied Party Supporter" & Patronage == "High", "Patron_Houthis", "Other"))
dat$Patron_Houthi <- as.factor(dat$Patron_Houthi)

# Rely job on state apparatuses
dat <- dat |> 
  mutate(Job_rely_state = ifelse(Q5_1 <= 5, "Rely job on state", "Others"))
dat$Job_rely_state <- as.factor(dat$Job_rely_state)

# Rely job on Tribal chief
dat <- dat |> 
  mutate(Job_rely_tribe = ifelse(Q5_1 == 6, "Rely job on tribe", "Others"))
dat$Job_rely_tribe <- as.factor(dat$Job_rely_tribe)

# Rely economic difficulty on state apparatuses
dat <- dat |> 
  mutate(economy_rely_state = ifelse(Q5_1 <= 5, "Rely economy on state", "Others"))
dat$economy_rely_state <- as.factor(dat$economy_rely_state)

# Rely security on state apparatuses
dat <- dat |> 
  mutate(Security_rely_state = ifelse(Q5_3 <= 5, "Rely security on state", "Others"))
dat$Security_rely_state <- as.factor(dat$Security_rely_state)

# Piety
dat <- dat |> 
  mutate(Piety = case_when(Q6 <= 4 ~ "Low",
                           Q6 >= 5 & Q6 <= 6 ~ "Middle",
                           Q6 >=7 ~ "High"))
dat$Piety <- as.factor(dat$Piety)

# Subjective conflict intensity
dat <- dat |> 
  mutate(Conflict_assessment = case_when(Q14 >= 4 ~ "Worse",
                                         Q14 == 3 ~ "Same",
                                         Q14 <= 2 ~ "Better"))
dat$Conflict_assessment <- as.factor(dat$Conflict_assessment)

# Subjective security perspective
dat <- dat |> 
  mutate(Future_expectation = case_when(Q15 >= 6 ~ "Worse",
                                        Q15 == 5 ~ "Same",
                                        Q15 <= 2 ~ "Better"))
dat$Future_expectation <- as.factor(dat$Future_expectation)

```

# whole model
```{r, fig.height = 5, dpi = 300}
f1 <- selected ~ Activity + Armament + Patron + Service
str(dat)

# MM
whole_mm <- cj(dat, f1, 
               id = id,
               estimate = "mm", 
               h0 = 0.5,
               alpha = 0.1)
whole_mm
plot(whole_mm, vline = 0.5, 
     size = 2)
write_html(whole_mm, file = "result/whole_mm_Yemen.html")

# amce
whole_amce <- cj(data = dat, f1,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 alpha = 0.1)
whole_amce
plot(whole_amce, vline = 0,
     size = 2)

whole_Yemen <- plot(whole_amce, vline = 0, size = 2) +
  ggplot2::labs(title = "Whole result of Yemen")
ggsave(filename = "result/whole_Yemen.png",
       plot = whole_Yemen,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

write_html(whole_amce, file = "result/whole_amce_Yemen.html")

# plot estimation difference
whole_mm$Estimate <- "MM"
whole_amce$Estimate <- "AMCE"
plot(rbind(whole_mm, whole_amce), feature_headers = FALSE, size = 2)+
  ggplot2::facet_wrap(~ Estimate, ncol = 2L)
```

# Ceasefire
```{r, fig.height = 5, dpi = 300}
# amce
ceasefire_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ceasefire,
                 alpha = 0.1)
ceasefire_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ceasefire,
                      alpha = 0.1)
plot(rbind(ceasefire_amce, ceasefire_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
ceasefire_mm <- cj(dat,
               selected ~ Activity + Armament + Patron + Service,
               id = id,
               estimate = "mm",
               by = ~ ceasefire,
               alpha = 0.1)
ceasefire_mm_diff <- cj(data = dat,
                    selected ~ Activity + Armament + Patron + Service,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ceasefire,
                    alpha = 0.1)
plot(rbind(ceasefire_mm, ceasefire_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(ceasefire_mm, group = "ceasefire", size = 2)

# plot
p1 <- plot(ceasefire_mm, group = "ceasefire", size = 2) +
  ggplot2::labs(title = "Yemen")
p2 <- plot(ceasefire_mm_diff, size = 2)+ 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
p3 <- ggarrange(p1,
                p2 + theme(axis.title.y = element_blank()),
                nrow = 1, ncol = 2, align = "hv",
                common.legend = FALSE,
                legend = "bottom")

ggsave(filename = "result/ceasfire_ym.png",
       plot = p3,
       width = 10,
       height = 4,
       units = "in", 
       dpi = 600)

# OLS
write_html(ceasefire_mm, file = "result/ceasefire_mm_ym.html")
write_html(ceasefire_mm_diff, file = "result/ceasefire_mm_diff_ym.html")

# anova
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ ceasefire)
```

# task
```{r, fig.height = 7, dpi = 300}
# amce
task_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ task)
task_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ task)
plot(task_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(rbind(task_amce, task_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# gender
```{r, fig.height = 5, dpi = 300}
# amce
gender_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Gender)
gender_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Gender)
plot(rbind(gender_amce, gender_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mmm
gender_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Gender)
gender_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Gender)
plot(rbind(gender_mm, gender_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# age group
```{r, fig.height = 7, dpi = 300}
# amce
Age_group_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Age_group)
Age_group_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Age_group)
plot(Age_group_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Age_group_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Age_group)
Age_group_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Age_group)
plot(Age_group_mm, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Age_group_mm, group = "Age_group", size = 2)
```

# Education
```{r, fig.height = 7, dpi = 300}
# amce
Education_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Education)
Education_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Education)
plot(Education_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Education_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Education)
Education_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Education)
plot(Education_mm, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Education_mm, group = "Education", size = 2)
```

# Income
```{r, fig.height = 7, dpi = 300}
# amce
Income_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Income)
Income_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Income)
plot(Income_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Income_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Income)
Income_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Income)
plot(Income_mm, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Sect
```{r, fig.height = 7, dpi = 300}
table(dat$Sect)
# amce
Sect_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Sect)
Sect_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Sect)
plot(rbind(Sect_amce, Sect_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Sect_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Sect)
Sect_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Sect)
plot(rbind(Sect_mm, Sect_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Sect_mm, group = "Sect", size = 2)
```

# Islamism
```{r, fig.height = 7, dpi = 300}
# amce
Islamism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Islamism)
Islamism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Islamism)
plot(rbind(Islamism_amce, Islamism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Islamism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Islamism)
Islamism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Islamism)
plot(rbind(Islamism_mm, Islamism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Islamism_mm, group = "Islamism", size = 2)

```

# Secularism
```{r, fig.height = 7, dpi = 300}
# amce
Secularism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Secularism)
Secularism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Secularism)
plot(rbind(Secularism_amce, Secularism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Secularism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Secularism)
Secularism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Secularism)
plot(rbind(Secularism_mm, Secularism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Secularism_mm, group = "Secularism", size = 2)

```

# Tribalism
```{r, fig.height = 7, dpi = 300}
# amce
Tribalism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Tribalism)
Tribalism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Tribalism)
plot(rbind(Tribalism_amce, Tribalism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Tribalism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Tribalism)
Tribalism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Tribalism)
plot(rbind(Tribalism_mm, Tribalism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Tribalism_mm, group = "Tribalism", size = 2)

```

# Arabism
```{r, fig.height = 7, dpi = 300}
# amce
Arabism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Arabism)
Arabism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Arabism)
plot(rbind(Arabism_amce, Arabism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Arabism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Arabism)
Arabism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Arabism)
plot(rbind(Arabism_mm, Arabism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward US
```{r, fig.height = 7, dpi = 300}
table(dat$US)
# amce
US_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ US)
US_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ US)
plot(rbind(US_amce, US_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
US_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ US)
US_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ US)
plot(rbind(US_mm, US_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(US_mm, group = "US", size = 2)
```

# Veiws toward Iran
```{r, fig.height = 5, dpi = 300}
table(dat$Iran)
# amce
Iran_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Iran,
                alpha = 0.1)
Iran_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Iran,
                     alpha = 0.1)
plot(rbind(Iran_amce, Iran_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Iran_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Iran,
              alpha = 0.1)
Iran_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Iran,
                   alpha = 0.1)
plot(rbind(Iran_mm, Iran_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# MM plot
plot(Iran_mm, group = "Iran", size = 2)
Iran_Yemen <- plot(Iran_mm, group = "Iran", size = 2) +
  ggplot2::labs(title = "Yemen")
ggsave(filename = "result/Iran_Yemen.png",
       plot = Iran_Yemen,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# MM difference plot
plot(Iran_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
Iran_diff_YM <- plot(Iran_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::labs(title = "Yemen")
ggsave(filename = "result/Iran_diff_YM.png",
       plot = Iran_diff_YM,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# anova
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ Iran)

# OLS result
write_html(Iran_mm, file = "result/Iran_mm_ym.html")
write_html(Iran_mm_diff, file = "result/Iran_mm_diff_ym.html")
```

# Veiws toward Saudi
```{r, fig.height = 7, dpi = 300}
table(dat$Saudi)
# amce
Saudi_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Saudi)
Saudi_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Saudi)
plot(rbind(Saudi_amce, Saudi_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Saudi_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Saudi)
Saudi_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Saudi)
plot(rbind(Saudi_mm, Saudi_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Turkey
```{r, fig.height = 7, dpi = 300}
table(dat$Turkey)
# amce
Turkey_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Turkey)
Turkey_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Turkey)
plot(rbind(Turkey_amce, Turkey_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Turkey_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Turkey)
Turkey_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Turkey)
plot(rbind(Turkey_mm, Turkey_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward UAE
```{r, fig.height = 7, dpi = 300}
table(dat$UAE)
# amce
UAE_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ UAE)
UAE_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ UAE)
plot(rbind(UAE_amce, UAE_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
UAE_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ UAE)
UAE_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ UAE)
plot(rbind(UAE_mm, UAE_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Israel
```{r, fig.height = 7, dpi = 300}
table(dat$Israel)
# amce
Israel_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Israel)
Israel_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Israel)
plot(rbind(Israel_amce, Israel_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Israel_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Israel)
Israel_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Israel)
plot(rbind(Israel_mm, Israel_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Axis of resistance
```{r, fig.height = 7, dpi = 300}
table(dat$Axis)
# amce
Axis_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Axis)
Axis_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Axis)
plot(rbind(Axis_amce, Axis_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Axis_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Axis)
Axis_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Axis)
plot(rbind(Axis_mm, Axis_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Axis_mm, group = "Axis", size = 2)
```

# Veiws toward UN
```{r, fig.height = 7, dpi = 300}
table(dat$UN)
# amce
UN_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ UN)
UN_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ UN)
plot(rbind(UN_amce, UN_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
UN_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ UN)
UN_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ UN)
plot(rbind(UN_mm, UN_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Houthi allied parties supporter
```{r, fig.height = 5, dpi = 300}
table(dat$Houthiallied_party_supporter)
# amce
Houthiallied_party_supporter_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Houthiallied_party_supporter,
                 alpha = 0.1)
Houthiallied_party_supporter_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Houthiallied_party_supporter,
                      alpha = 0.1)
plot(rbind(Houthiallied_party_supporter_amce, Houthiallied_party_supporter_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Houthiallied_party_supporter_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Houthiallied_party_supporter,
                 alpha = 0.1)
Houthiallied_party_supporter_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Houthiallied_party_supporter,
                      alpha = 0.1)
plot(rbind(Houthiallied_party_supporter_mm, Houthiallied_party_supporter_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Houthiallied_party_supporter_mm, group = "Houthiallied_party_supporter", size = 2)
```

# Patronage
```{r, fig.height = 5, dpi = 300}
# amce
Patronage_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Patronage,
                 alpha = 0.1)
Patronage_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Patronage,
                      alpha = 0.1)
plot(rbind(Patronage_amce, Patronage_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Patronage_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Patronage,
                 alpha = 0.1)
Patronage_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Patronage,
                      alpha = 0.1)
plot(rbind(Patronage_mm, Patronage_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Patronage_mm, group = "Patronage", size = 2)
```

# High patronage and Houthi allied supporter
```{r, fig.height = 5, dpi = 300}
# amce
Patron_Houthi_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Patron_Houthi,
                 alpha = 0.1)
Patron_Houthi_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Patron_Houthi,
                      alpha = 0.1)
plot(rbind(Patron_Houthi_amce, Patron_Houthi_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Patron_Houthi_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Patron_Houthi,
                 alpha = 0.1)
Patron_Houthi_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Patron_Houthi,
                      alpha = 0.1)
plot(rbind(Patron_Houthi_mm, Patron_Houthi_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Patron_Houthi_mm, group = "Patron_Houthi", size = 2)
```

# Rely job on state apparatuses
```{r, fig.height = 5, dpi = 300}
table(dat$Job_rely_state)
# amce
Job_rely_state_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Job_rely_state)
Job_rely_state_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Job_rely_state)
plot(rbind(Job_rely_state_amce, Job_rely_state_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Job_rely_state_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Job_rely_state)
Job_rely_state_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Job_rely_state)
plot(rbind(Job_rely_state_mm, Job_rely_state_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Rely job on Tribal chief
```{r, fig.height = 5, dpi = 300}
table(dat$Job_rely_tribe)
# amce
Job_rely_tribe_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Job_rely_tribe)
Job_rely_tribe_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Job_rely_tribe)
plot(rbind(Job_rely_tribe_amce, Job_rely_tribe_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Job_rely_tribe_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Job_rely_tribe)
Job_rely_tribe_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Job_rely_tribe)
plot(rbind(Job_rely_tribe_mm, Job_rely_tribe_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Rely economic dificulty on state apparatuses
```{r, fig.height = 5, dpi = 300}
table(dat$economy_rely_state)
# amce
economy_rely_state_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ economy_rely_state)
economy_rely_state_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ economy_rely_state)
plot(rbind(economy_rely_state_amce, economy_rely_state_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
economy_rely_state_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ economy_rely_state)
economy_rely_state_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ economy_rely_state)
plot(rbind(economy_rely_state_mm, economy_rely_state_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Rely security on state apparatuses
```{r, fig.height = 5, dpi = 300}
table(dat$Security_rely_state)
# amce
Security_rely_state_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Security_rely_state)
Security_rely_state_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Security_rely_state)
plot(rbind(Security_rely_state_amce, Security_rely_state_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Security_rely_state_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Security_rely_state)
Security_rely_state_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Security_rely_state)
plot(rbind(Security_rely_state_mm, Security_rely_state_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Piety
```{r, fig.height = 7, dpi = 300}
table(dat$Piety)
# amce
Piety_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Piety)
Piety_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Piety)
plot(rbind(Piety_amce, Piety_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Piety_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Piety)
Piety_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Piety)
plot(rbind(Piety_mm, Piety_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Piety_mm, group = "Piety", size = 2)
```

# Conflict_assessment
```{r, fig.height = 7, dpi = 300}
table(dat$Conflict_assessment)
# amce
Conflict_assessment_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Conflict_assessment)
Conflict_assessment_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Conflict_assessment)
plot(rbind(Conflict_assessment_amce, Conflict_assessment_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Conflict_assessment_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Conflict_assessment)
Conflict_assessment_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Conflict_assessment)
plot(rbind(Conflict_assessment_mm, Conflict_assessment_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Conflict_assessment_mm, group = "Conflict_assessment", size = 2)
```

# Future_expectation
```{r, fig.height = 7, dpi = 300}
table(dat$Future_expectation)
# amce
Future_expectation_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Future_expectation)
Future_expectation_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Future_expectation)
plot(rbind(Future_expectation_amce, Future_expectation_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Future_expectation_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Future_expectation)
Future_expectation_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Future_expectation)
plot(rbind(Future_expectation_mm, Future_expectation_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Future_expectation_mm, group = "Future_expectation", size = 2)
```

# respondant location
```{r}
table(dat$LocationLatitude, dat$LocationLongitude)
```

# frequency and propotion
```{r}
plot(cj_freqs(dat,
              f1,
              id = ~ respondentIndex))
```


# interaction: Ceasefire*Houthiallied_party_supporter
```{r, fig.height = 5, dpi = 300}
dat$Houthis_ceasefire <- interaction(dat$Houthiallied_party_supporter, dat$ceasefire, sep = "_")
Houthi_ceasefire_mm <- cj(dat,
               selected ~ Houthis_ceasefire,
               id = id,
               estimate = "mm",
               alpha = 0.1)
plot(Houthi_ceasefire_mm, vline = 0.5, size = 2)

# plot
Houthi_ceasefire <- plot(Houthi_ceasefire_mm, vline = 0, size = 2)+
  labs(title = "Yemen")
ggsave(filename = "result/Houthi_ceasefire.png",
       plot = Houthi_ceasefire,
       width = 5,
       height = 3,
       units = "in", 
       dpi = 600)

# OLS result
write_html(Houthi_ceasefire_mm, file = "result/Houthi_ceasefire_mm.html")
```


# interaction
```{r}
table(dat$Iran, dat$Houthiallied_party_supporter)
cor(as.numeric(dat$Iran), as.numeric(dat$Houthiallied_party_supporter))

# subset ceasefire and shia party supporter
pre_supporter <- dat |>
  filter(ceasefire == "Pre")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ Houthiallied_party_supporter,
     alpha = 0.1)
plot(pre_supporter, group = "Houthiallied_party_supporter", 
     size = 2,
     vline = 0.5)

post_supporter <- dat |>
  filter(ceasefire == "Post")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ Houthiallied_party_supporter,
     alpha = 0.1)
plot(post_supporter, group = "Houthiallied_party_supporter", 
     size = 2,
     vline = 0.5)

plot(rbind(pre_supporter, post_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)



pre_supporter_diff <- dat |>
  filter(ceasefire == "Pre")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ Houthiallied_party_supporter,
     alpha = 0.1)
plot(pre_supporter_diff, group = "Houthiallied_party_supporter", 
     size = 2,
     vline = 0.5)

post_supporter_diff <- dat |>
  filter(ceasefire == "Post")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ Houthiallied_party_supporter,
     alpha = 0.1)
plot(post_supporter_diff, group = "Houthiallied_party_supporter", 
     size = 2,
     vline = 0.5)

plot(rbind(pre_supporter, pre_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(rbind(post_supporter, post_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


plot(rbind(pre_supporter_diff, post_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

###
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
supporter <- dat |>
  filter(Houthiallied_party_supporter == "Houthis Allied Party Supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)

non_supporter <- dat |>
  filter(Houthiallied_party_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)
plot(rbind(supporter, non_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

supporter_diff <- dat |>
  filter(Houthiallied_party_supporter == "Houthis Allied Party Supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter_diff, group = "ceasefire", 
     size = 2)

non_supporter_diff <- dat |>
  filter(Houthiallied_party_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter_diff, group = "ceasefire", 
     size = 2)


# interaction
supporter$Supporter <- "Supporter"
non_supporter$Supporter <- "Non_supporter"
interaction <- rbind(supporter, non_supporter)

H2_app_ym <- plot(interaction, group = "Supporter", size = 2) +
  ggplot2::facet_wrap(~ BY, ncol = 2) +
  labs(title = "Yemen",
       shape = "Hamas Supporter",
       labels = c("Non-suporter", "Supporter"))+
  scale_color_hue(labels = c("Non-suporter", "Supporter"))
H2_app_ym

ggsave(filename = "result/H2_app_ym.png",
       plot = H2_app_ym,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)  

# OLS result
write_html(interaction, file = "result/H2_app_ym.html")

# diff
supporter_diff$Supporter <- "Supporter"
non_supporter_diff$Supporter <- "Non_supporter"
interaction_diff <- rbind(supporter_diff, non_supporter_diff)

H1_app_ym <- plot(interaction_diff, group = "Supporter", size = 2) +
  ggplot2::facet_wrap(~ BY, ncol = 2) +
  labs(title = "Yemen",
       labels = c("Non-suporter", "Supporter"))+
  scale_color_hue(labels = c("Non-suporter", "Supporter"))
H1_app_ym

ggsave(filename = "result/H1_app_ym.png",
       plot = H1_app_ym,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)  

```
